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Abstract 


Sub-mesoscale  structures  such  as  oceanic  fronts  can  play  a  very  important  role  in  the 
transfer  of  properties,  tracers,  momentum  and  energy  between  the  surface  ocean  and  the 
upper  thermocline.  The  ability  to  correctly  represent  the  fronts  in  ocean  circulation  models  is 
important  for  mixed  layer  and  thennocline  simulations  as  well  as  large  scale  circulations. 

Numerical  experiments  are  conducted  using  hydrostatic  (HY)  and  nonhydrostatic 
(NH)  models  to  address  the  relevance  of  NH  effects  on  the  evolution  of  density  fronts  and 
the  development  of  meso-  and  submeso-scale  vertical  motions.  Model  results  indicate  that 
even  though  the  NH  simulations  give  stronger  upwelling/downwelling  than  the  HY 
simulations,  the  characteristics  of  submesoscale  structures  are  very  similar  between  them, 
suggesting  that  the  HY  model  are  capable  to  produce  correct  dynamical  response  for  oceanic 
fronts. 

Model  resolution  is  found  to  be  critical  for  frontal  simulations  for  both  hydrostatic 
and  nonhydrostatic  models,  and  100  meter  resolution  is  adequate  for  frontal  simulations. 


v 


1.  Introduction 


Oceanic  submesoscale  structures  usually  have  a  horizontal  scale  of  1  to  10  km.  They 
are  inevitably  crucial  to  bridging  the  meso  and  smaller  scales  (Thomas  et.  al.  2008).  At  the 
submesoscale,  the  nonhydrostatic  process  start  to  kick  in,  but  may  or  may  not  be  significant 
depending  on  what  phenomenon  we  are  looking  at.  Internal  tides  are  examples  where 
nonhydrostatic  processes  have  significant  influence,  while  the  oceanic  fronts  can  still  be 
considered  hydrostatic  given  their  horizontal  and  vertical  scales  (Mahadevan  2006). 

Fronts  are  very  common  to  the  ocean  (Ullman  and  Cornillon  1999),  and  associated 
with  narrow  geostrophic  jets  and  strong  shear  (Ferrari  and  Rudnick,  2000).  Sub-mesoscale 
structures  such  as  fronts  are  thought  to  be  instrumental  in  transferring  energy  and  properties 
from  the  largely  adiabatic  mesoscale  flow  field  to  a  scale  where  mixing  occurs  (McWilliams 
2003,  Molemarker  et  al.  2005,  Thomas  et.  al.  2008).  They  can  facilitate  the  exchange  of 
properties  between  the  surface  ocean  and  upper  thennocline.  The  intensity  of  the  frontal,  up 
and  down-welling  motions  will  set  the  rates  of  exchange  in  the  vertical,  particularly  beyond 
the  reach  of  the  mixed  layer.  The  ability  to  correctly  represent  them  in  the  ocean  models  is 
important  for  mixed  layer  and  thennocline  simulations  as  well  as  large  scale  circulations. 
The  location  of  density  fronts  and  the  depth  of  mixed  layer  and  thennocline  are  also  very 
important  for  acoustic  perfonnance  predictions. 

Marshall  et  al  (1997)  argue  that  the  hydrostatic  models  presumably  begin  to  break 
down  somewhere  between  10  and  1  km  resolution,  as  the  resolved  horizontal  scale  of  the 
motion  becomes  comparable  with  its  vertical  scale.  However,  nonhydrostatic  models  have 
not  shown  much  advantage  on  simulating  oceanic  fronts  vs.  hydrostatic  models  at 
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resolutions  from  0.25  km  to  1  km.  Mahadevan  (2006)  argues  that  greater  resolution  is 
required  to  enter  the  regime  where  the  differences  become  apparent.  Although  it  is  well 
established  that  higher  resolution  can  produce  more  realistic  modeling  of  small-scale 
features,  the  importance  of  grid  resolution  in  the  vicinity  of  strong  density  fronts  has  not  yet 
been  investigated.  In  this  study,  hydrostatic  and  nonhydrostatic  models  are  used  to  simulate 
oceanic  salinity  fronts  at  resolutions  range  from  50  meters  to  1km.  The  model  details  are 
given  in  section  2;  section  3  outlines  the  design  of  numerical  experiments  conducted;  in 
section4,  the  characteristics  of  the  submesoscale  motions  are  assessed  in  these  simulations 
and  the  underlying  mechanisms  and  conditions  for  the  existence  and  evolution  of  the 
submesoscale,  frontal  processes  and  the  vertical  circulation  associated  with  them  are 
investigated;  the  conclusion  is  given  in  section  5. 


2,  Method 


2.1  The  Nonhydro  static  model 

NHWAVE  is  a  non-hydrostatic  wave  and  circulation  model  developed  at  University 


of  Delaware  (Ma  et  al.  2012,  Shi  et  ah,  2015).  Here,  we  only  list  the  equations  related  to  the 


present  study.  The  incompressible  Navier-Stokes  equations  can  be  written  as 
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where  u(x,  y,  z,  t),  v(x,  y,  z,  t),  w(x,  y,  z,  t)  are  flow  velocity  components,  and  v  =  (u,  v,  w); 
p  is  water  density;  p(x,  y,  z,  t)  is  pressure;  gx,  gy,  gz  are  gravitational  force  components;  x,j, 
(i,  j)  =  x,  y,  z  are  turbulent  stress  components  which  can  be  modeled  by  a  turbulence  closure 
model.  In  the  present  study,  the  standard  Smagorinsky  LES  model  is  used.  The  conservation 
equation  for  a  passive  tracer,  for  example,  salinity  S,  can  be  written  as 


dS  dS  dS  dS  _ 

—  +  U  —  +  v  —  +  w—  =  ksV2S 

at  ox  dy  dz 


(5) 


where  ks  is  the  diffusivity  of  salinity.  The  fluid  density  p  is  calculated  from  salinity  and 
temperature  based  on  the  Equation  of  State  of  seawater.  However,  a  constant  temperature  is 
used  in  the  following  test  cases.  The  Boussinesq  approximation  is  employed  for  fonnulating 
the  baroclinic  pressure  in  the  momentum  equations. 

NHWAVE  uses  a  combined  finite-volume  and  finite-difference  scheme  with  a 


Godunov-type  method  for  the  spatial  discretization  and  a  two-stage  (second-order)  Strong 
Stability-Preserving  (SSP)  Runge-Kutta  scheme  for  the  time  integration.  A  two-step 
predictor-corrector  projection  method  is  applied  to  solve  the  dynamic  pressure.  The  high 
performance  HYPRE  software  library  is  used  to  solve  the  pressure  Poisson  equation.  The 
model  is  parallelized  with  Message  Passing  Interface  (MPI). 


2.2  Modification  to  NHWAVE 


NHWAVE  is  designed  for  coastal  simulations  with  small  domains.  The  earth 
rotation  effect  is  not  considered.  To  use  NHWAVE  for  ocean  front  simulations,  Coriolis 
effect  is  added  to  the  model  momentum  equations: 
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Where  /  =  2fl  sin  <fi  and  b  =  2  ficos  <p  ((/>  is  the  latitude)  are  the  components  of  Coriolis 
acceleration  in  the  horizontal  (x,  y)  and  vertical  (z)  directions.  Note  that  the  Coriolis  terms 
are  based  on  the  Cartesian  coordinates  with  y  pointing  to  north. 


2.3  Hydrostatic  model  vs.  Nonhydrostatic  model 

In  the  nonhydrostatic  model,  all  terms  in  the  incompressible  Navier  Stokes  equations 
(4)  -  (8)  are  retained,  while  in  the  hydrostatic  model,  all  the  tenns  underlined  by  the  curly 
brackets  are  neglected.  The  vertical  velocity  w  is  diagnosed  from  the  continuity  equation  (4) 
in  the  hydrostatic  model  rather  than  obtained  prognostically  from  equation  (8)  as  the 
nonhydrostatic  model  does. 

By  reducing  the  vertical  momentum  equation  to  a  statement  of  hydrostatic  balance, 
the  hydrostatic  (HY)  model  differs  from  the  nonhydrostatic  (NH)  model  due  to  the  exclusion 
of  the  nonhydro  static  pressure  gradient,  and  the  Coriolis  acceleration  terms  that  contain  the 
component  of  the  earth’s  angular  velocity,  2f2cosc|),  in  a  direction  tangential  to  the  surface  of 
the  earth.  NHWAVE  is  coded  to  easily  switch  between  HY  and  NH  mode  through  input  file 
specifications. 
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2.4  Geostrophic  balance 

Both  HY  and  NH  model  experiments  are  initialized  using  geostrophic  balance.  If  we 
separate  the  pressure  p  into  its  HY  component  phy  and  NH  component  pnh,  the  geostrophic 
balance  for  the  NH  model  can  be  written  as 


f^  =  dJ^+yd-T~Robw 


-f^  =  dJ^  +  r-e f 

Y  dPNH 


bu 


(10) 

(11) 

(12) 


y  S  dz 

where,  y  is  the  ratio  of  NH  to  HY  pressure  variations,  5  is  the  depth  to  length  aspect  ratio, 
and  R0  =  U/QL  is  the  Rossby  number.  In  the  NH  model,  b  is  non-zero  and  y  =  §.  The  HY 
case  is  obtained  by  setting  y  =  0  and  b  =  0. 


3.  Numerical  Experiments 

NHWAVE  is  run  in  both  HY  and  NH  modes  in  four  sets  of  experiments  with  the 
resolution  of  500,  250,  100,  and  50  meters  respectively. 

In  all  experiments,  the  model  domain  is  set  up  to  be  16  km  in  the  east- west  direction, 
24  km  in  the  north-south  direction  and  800  meters  in  the  vertical  direction.  An  idealized 
front  is  initialized  using  geostrophic  balance  as  shown  in  Figure  1.  Periodic  boundary 
conditions  are  applied  in  the  x  direction. 

The  domain  is  initialized  with  a  sharp  lateral  south-to-north  (S-N)  density  gradient 
based  on  the  variation  of  salinity  in  the  upper  layers,  such  that  the  southern  half  of  the 
domain  has  lighter  water.  This  front  is  representative  of  the  deep,  semi-permanent  fronts  that 
are  observed  in  the  ocean  (Rudnick,  1996)  and  extends  to  a  depth  of  250  meters.  The 
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vertical  stratification  is  prescribed  based  on  observed  open  ocean  profile  in  the  North 
Atlantic.  It  overlies  a  weakly  stratified  region  and  a  nearly  homogeneous  deeper  layer. 
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Figure  1.  Across-channel  cross  section  showing  the  sea  surface  elevation  (top),  potential 
density  (middle)  and  along  channel  velocity  in  the  upper  400  m  that  were  used  to  initialize 
the  model  with  an  upper  ocean  front. 
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A  S-N  salinity  variation  is  then  superimposed  on  this  vertical  stratification.  In  the  upper 
50m  AS  =  0.2psu;  between  50m  and  250m  AS  =  0.2(250-z)/200,  where  z  is  the  depth  in  m; 
and  below  250m,  there  is  no  horizontal  density  variation.  The  salinity  variation  in  the  S-N 
direction  is  distributed  according  to  (±AS/2)(l-exp(-yc/2))/(l+exp(-yc/2)),  where  yc  is  the 
distance  in  km  from  the  channel  center  (northward  is  positive),  and  the  ±  sign  are  used 
according  as  yc  is  positive  or  negative.  The  sea  surface  elevation  is  varied  correspondingly 
over  the  same  frontal  region  by  3.8cm,  being  higher  in  the  southern  region.  Associated  with 
the  S-N  density  front  is  a  west-to-east  (W-E)  geostrophic  jet.  The  model  velocities  (and  the 
NH  pressure  in  the  NH  model)  are  initialized  to  satisfy  the  model’s  respective  NH  or  HY 
balanced  state  as  described  in  section  2.3.  The  model  is  allowed  to  evolve  from  this  state 
into  a  meandering  front,  generating  smaller  scale  features,  when  perturbed  by  weak  winds. 
A  constant  westerly  wind  stress  of  0.2Nm'2  (corresponding  to  lOm/s  winds)  is  applied  to  the 
entire  model  domain  continuously  for  50  days. 

4.  Results 

In  this  section,  we  analyze  the  numerical  results  to  address  the  relevance  of  NH 
effects  on  the  evolution  of  density  fronts  and  the  development  of  meso-  and  submeso-scale 
vertical  motions.  The  numerical  resolution  required  for  modeling  the  vertical  velocity  at  the 
fronts  is  also  investigated.  Based  on  these  results,  we  will  assess  whether  hydrostatic  models 
are  capable  on  simulating  these  features. 
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4.1  HY  simulation  with  500  m  resolution 


First,  the  evolution  of  the  front  is  examined  for  the  hydrostatic  model  with  500  m 
resolution.  As  the  front  evolves,  the  surface  density  field  develops  meanders.  On  day  31, 
submesoscale,  vertical  circulation  cells  are  formed  close  to  the  surface  along  the  boundary 
of  the  meandering  front  (Figure  2,  upper  panels).  The  vertical  velocity  field  at  both  20  and 
70  meter  depth  show  submesoscale  features  with  the  upwelling/downwelling  velocities 
being  stronger  at  70  m  depth  where  the  base  of  the  mixed  layer  is. 
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Figure  2.  The  surface  density  and  the  vertical  velocity  field  (m/s)  shown  at  20  m  and  70  m 
depth  on  day  31  from  the  HY  (upper  panels)  and  NH  (lower  panels)  model  forced  with 
constant  westerly  wind  stress  of  0.2  N/m2. 
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Stronger  perturbation  to  the  front  is  observed  for  the  NH  simulation.  From  the 
surface  density  plot,  we  can  clearly  see  an  eddy  being  formed  and  pinched  off  the 
meandering  front.  The  magnitude  of  the  upwelling/downwelling  velocities  have  the  similar 
magnitude  as  the  HY  simulation,  but  the  submesoscale  features  in  the  vertical  velocity  field 
are  elongated  in  the  east-west  direction.  This  may  due  to  the  fact  that  the  growth  rate  for 
baroclinic  instability  in  the  N  H  and  HY  cases  are  different  (M  olemakeret.  al.  2005). 


x  (km) 


x  (km) 


x  (km) 


Figure  3.  The  surface  density  and  the  vertical  velocity  field  (m/s)  shown  at  20  m  and  70  m 
depth  on  day  50  from  the  HY  (upper  panels)  and  NH  (lower  panels)  model  forced  with 
constant  westerly  wind  stress  of  0.2  N/m2. 
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On  day  50,  the  density  and  vertical  velocity  fields  become  very  similar  between  the 
HY  and  NH  simulations  (Figure  3).  The  density  gradients  of  the  front  are  much  smaller,  and 
the  upwelling/downwelling  velocities  are  much  weaker.  The  similarity  between  the  HY  and 
NH  simulations  suggests  that  these  submesoscale  features  developed  along  the  fronts  are 
HY  phenomenon  or  dominated  by  HY  dynamics. 

4.2  Comparison  of  HY  andNH  solutions  at  varies  grid  resolutions. 

On  account  of  the  nonlinear  evolution  of  the  flow,  the  velocities  and  density 
distributions  are  different  between  the  HY  and  NH  simulations  at  500  m  resolution. 
However,  it  is  difficult  to  identify  characteristic  differences  between  the  NH  and  HY 
solutions.  The  differences  between  NH  and  HY  models  are  more  significant  at  the  onset  of 
baroclinic  instability  because  the  evolution  of  the  instability  differs  between  the  two  models 
(Mahadevan  2006).  The  wavy  defonnations  of  surface  isopycnals  tends  to  be  more  angularly 
oriented  in  the  direction  of  the  mean  zonal  flow  in  the  NH  case.  But  with  time,  the  HY  and 
NH  solutions  appear  statistically  similar.  Thus,  in  Figure  4,  frontal  features  from  the  NH  and 
HY  models  at  various  horizontal  resolution  (500  m,  250  m,  and  100  m)  are  compared  at  an 
earlier  stage  of  the  simulation  (5  days). 
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Figure  4.  The  near  surface  vertical  velocity  at  the  upper  ocean  front  5  days  after  model 
initialization  in  the  NY  and  NH  models  at  grid  resolution  of  500,  250,  and  100  m.  The  color 
bar  show  the  magnitude  of  w  in  m/s.  Note  that  the  range  of  color  bar  for  (c)  and  (f)  are  larger 
due  to  the  stronger  upwel I i ng/dow nwel I i ng  in  the  100m  resolution  simulations. 


Instability  seems  to  develop  much  slower  at  lower  resolution  (500  m).  Only  weak 
meanders  along  the  front  are  observed  at  day  5  for  both  the  HY  and  NH  simulations.  The 
structure  and  magnitude  of  the  near  surface  vertical  velocity  are  very  similar  between  them. 
When  the  resolution  is  doubled  (250  m),  submesoscale  structures  start  to  develop.  Large 
areas  of  upwel  ling  are  accompanied  by  long  strips  of  stronger  dow  nwel  ling  along  the  edges 
of  the  meandering  front  where  outcropping  surface  isopycnals  are.  The 
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upwelling/downwelling  velocities  are  at  the  similar  magnitude  for  the  HY  and  NH 
simulations.  Although  the  strong  up  welling  regions  appear  at  different  locations  due  to  the 
random  nature  of  the  instability,  the  results  are  not  considered  statistically  different  between 
the  HY  and  NH  simulations. 

As  we  increase  the  model  resolution  to  100  m,  more  submesoscale  structures  are 
observed.  The  width  of  the  downwelling  zone  is  narrower  and  sharper,  and  their  magnitude 
increased  by  more  than  30%  compare  to  the  250  m  resolution  simulations.  The  NH 
simulation  gives  stronger  upwelling/downwelling  velocities  than  the  HY  simulation,  but  the 
HY  simulation  can  still  produce  the  similar  dynamical  structure  as  the  NH  simulation. 
Further  increase  the  resolution  to  50  m  does  not  alter  the  results  much  (not  shown), 
suggesting  that  the  frontal  submesoscale  features  can  be  adequately  resolved  at  100  m 
resolution. 

The  domain  wide  maximum  surface  upwelling/downwelling  velocities  from  the  NH 
and  HY  model  runs  at  250m  and  100m  resolution  are  compared  in  Figure  5.  The  region 
considered  excludes  a  strip  along  the  solid  northern  and  southern  boundaries  of  the  domain, 
where  upwelling  is  seen  due  to  boundary  effects. 
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Maximum  surface  upwelling  velocity 


Figure  5.  The  maximum  surface  upwel I i ng/downwel I i ng  velocity  plotted  against  time  for  the 
250m  and  100m  resolution  simulations  for  both  HY  and  NH  models.  Downwelling  is  more 
intense  than  upwelling,  and  increasing  grid  resolution  increases  the  most  intense  downward 
velocities. 
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Maximum  upwelling  velocity  at  70m  depth 


Figure  6.  The  maximum  upwel I i ng/downwel I i ng  velocity  at  70m  depth  plotted  against  time 
for  the  250m  and  100m  resolution  simulations  for  both  HY  andNH  models. 


N  H  models  start  to  develop  instability  earlier  than  HY  models,  but  the  magnitude  of 
the  disturbance  eventually  converge  between  the  HY  and  NH  models  in  time.  Thus,  after 
many  days  of  model  simulations,  it  is  difficult  to  find  systematic  differences  in  the  vertical 
velocities  between  the  NH  and  HY  models.  Also  notice  that  the  variation  of  downwelling 
maxima  at  the  mixed  layer  base  show  big  difference  between  HY  and  NH  simulations  for 
both  resolutions. 
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At  the  ocean  surface,  the  downwelling  maximum  become  more  intense  with 
increasing  grid  resolution,  but  the  upwelling  maxima  are  not  affected  by  the  model 
resolution  much  except  during  the  first  20  hours  or  so.  While  at  the  base  of  the  mixed  layer 
(~70m  depth),  both  upwelling  and  downwelling  maxima  are  affected  by  model  resolution. 
The  higher  resolution  simulations  give  stronger  upwelling/downwelling  (Figure  6). 

5.  Discussion  and  Conclusion 

Sub-mesoscale  structures  such  as  fronts  can  play  a  very  important  role  in  the  transfer 
of  properties,  tracers,  momentum  and  energy  between  the  surface  ocean  and  the  upper 
thermocline.  They  are  thought  to  be  crucial  in  bridging  the  meso  and  smaller  scales  (Thomas 
et.  al.  2008).  The  ability  to  correctly  represent  the  fronts  in  ocean  circulation  models  is 
important  for  mixed  layer  and  thennocline  simulations  as  well  as  large  scale  circulations. 
Since  they  are  not  fully  three-dimensional  and  nonhydrostatic,  it  is  important  to  address  the 
question  whether  they  can  be  adequately  simulated  in  hydrostatic  models  due  to  the  fact  that 
hydrostatic  ocean  circulation  models  are  still  the  most  widely  used  tool  in  the  ocean 
community. 

In  this  study,  numerical  experiments  are  conduced  using  hydrostatic  (HY)  and 
nonhydrostatic  (NH)  models  on  the  response  of  oceanic  fronts  to  surface  wind  forcing  to 
provide  a  descriptive  picture  of  submesoscale,  frontal  processes  and  the  vertical  circulation 
associated  with  them.  Model  results  indicate  that  even  though  the  NH  simulations  give 
stronger  upwelling/downwelling  than  the  HY  simulations,  the  characteristics  of 
submesoscale  structures  are  very  similar  between  them,  suggesting  that  the  HY  model  are 
capable  to  produce  correct  dynamical  response  for  oceanic  fronts.  Model  resolution  is  found 
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to  be  critical  for  frontal  simulations  for  both  hydrostatic  and  nonhydrostatic  models,  and  100 
meter  resolution  is  adequate  for  frontal  simulations. 

Further  study  of  the  frontal  processes  is  limited  by  the  MPI  design  of  NHWAVE. 
Although  it  supports  multi  CPU  computations,  the  model  grid  is  only  decomposed  in  the 
horizontal.  Thus  the  number  of  CPUs  can  be  used  are  constrained  by  the  number  of 
horizontal  grid  points.  Our  tests  suggest  that  one  CPU  need  to  have  about  100  horizontal 
grid  points  in  order  to  optimized  the  model  speed  for  the  chosen  model  grid.  Further 
increase  CPU  numbers  will  actually  slow  down  the  model  due  to  increased  communications 
among  CPUs.  Also,  the  scale  of  this  model  is  poor.  Increase  model  domain  size  by  4  times 
will  slow  down  the  model  by  more  than  10  times.  For  example,  the  NH  experiment  with 
100m  resolution  can  only  compute  ~6  physical  hours  in  1  day.  These  limitation  has 
tremendously  slowed  down  our  experiments  and  forced  us  to  use  very  small  model  domains. 
As  a  result,  we  can  only  resolve  about  half  a  wave  length  across  the  front,  which  makes  it 
impossible  for  us  to  further  analyze  the  mechanisms  underlie  the  formation  of  the 
mesolscale  motions,  for  example:  what  set  the  length  scale  of  the  mesoscale  structures  and 
affect  their  distribution  in  space  and  time,  and  what  affect  the  intensity  of  the  vertical 
velocity,  etc. 
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